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The existence of a slowly always-decreasing probability density for the recurrence times of earth- 
quakes implies that the occurrence of an event at a given instant becomes more unlikely as time 
since the previous event increases. Consequently, the expected waiting time to the next earthquake 
increases with the elapsed time, that is, the event moves away to the future fast. We have found 
direct empirical evidence of this counterintuitive behavior in two worldwide catalogs as well as in 
diverse local catalogs. Furthermore, the phenomenon can be well described by universal scaling 
functions. 

PACS numbers: 91.30.Dk, 05.65.+b, 89. 75. Da, 64.60.Ht 



Many probability distributions have been proposed to 
account for the recurrence time of earthquakes 
which is the time interval between successive earthquakes 
in a certain region. When aftershocks and mainshocks 
are considered together as a part of essentially one unique 
process J5|,[| , we have determined that a universal scaling 
law describes the probability density D(t) of the recur- 
rence time. In this way, for events of magnitude M above 
a certain threshold M c in a given spatial area (whose lim- 
its do not need to depend on the tectonic background), 
D(t) scales with the rate of seismic activity R in the area 
as 

D{t) = Rf(Rr), 

where R is defined as the mean number of earthquakes 
(with M > M c ) per unit time and / is a universal scal- 
ing function. This scaling is in fact the hallmark of the 
self-similarity of seismicity in the space-time-magnitudc 
domain. Among several possible scaling functions, the 
best fit is obtained from a (truncated) gamma distribu- 
tion, 
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with 7 the shape parameter, a the scale parameter, C 
a correction to normalization (due to the truncation of 
the distribution close to zero), and IX7) the gamma func- 
tion. When 7 < 1, / turns out to be a decreasing power 
law, accelerated by an exponential factor in the long-time 
limit. 

We analyze many regions and M c — values from two 
worldwide catalogs (the NEIC-PDE and the Significant 
Earthquake Database from NOAA at NEIC 0) and from 
several regional catalogs (Southern California, Japan, 
New Zealand, New Madrid (USA), the Iberian Penin- 
sula, and the British Islands ||). Each analyzed region is 
delimited by two meridians and two parallels || , with lin- 
ear size spanning from 0.16° (about 18 km) to the whole 
globe (~ 20 • 10 3 km), covering a large variety of tectonic 



environments, whereas the considered magnitudes range 
from larger than 1.5 to larger than 7.5 (this is about a 
factor 10 9 in the minimum released energy). 

Except for a 360° x 180°— region, which covers the en- 
tire globe, the rest of the regions are defined by a window 
of L degrees in longitude and L degrees in latitude. The 
coordinates (x,y) of the west-south corner of these re- 
gions can be obtained from the vector (k x ,k y ) at Fig. 
l's labels as x — x min + k x L, y = y rnin + k y L, with 
(x mm ,y mm ) = (-180°, -90°), (-123°, 30°), (127°, 27°), 
(160°, -60°), (-91°, 35°), (-20°, 30°), and (-10°, 45°) 
for the worldwide, Southern California, Japan, New 
Zealand, New Madrid, Iberian Peninsula, and British 
Island catalogs, respectively. The periods of study are 
(in years A.D. including the last one) 1973-2002 for the 
NEIC catalog, 1897-1970 for the NOAA one, 1988-1991, 
1995-1998, and 1984-2001 for Southern California (de- 
noted as SC88, SC95, and SC84), and 1995-1998, 1996- 
2001, 1975-2002, 1993-1997, and 1991-2001, for the rest 
of catalogs (in the same order as before). The regions and 
times of observation are selected in order that a period 
of stationary seismic activity is included, this means that 
aftershock sequences should not have too much weight in 
the seismicity of the region. When this is not the case 
(i.e., for very large aftershock activity) our analysis is 
still valid, but the scaling with the mean rate has to be 
replaced by a scaling with the instantaneous rate. 

A maximum-likelihood estimation of the parameters 
using the rescaled recurrence time Rt for all the regions 
and M c 's studied gives 7 = 0.74 ± 0.05, and a — 1.23 ± 
0.15, which yields a coefficient of variation CV ~ 1.2. 
The constant C is determined from the normalization 
condition given the minimum value for which the gamma 
distribution holds; for 9 > 0.05, C = 1.10 ±0.10 (see Fig. 
1 (a)). 

The results of the fit are shown in Fig. 1 (a) using the 
survivor function, which is defined as S(t) = Prober' > 
t] = f°° D{r')dr' (where r' is a generic label for the 
recurrence time, while r refers to a particular value of 
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the same quantity). It is straightforward to obtain that, 
in our case, S(t) should also verify a scaling relation, 
S(t) = G(Rt), with G{9) = CQ^O/a), and Q 7 (6»/a) the 
complement of the incomplete gamma function @,|nj. 
The total agreement between these equations and the 
measured distributions is clear from the data collapse 
and the fitting curve in the plot, for intermediate and 
long recurrence times. The accuracy of the scaling law 
and the gamma fit is guaranteed as the seismic activity 
is stationary in this range of recurrence times. On the 
contrary, short times are usually not free of disturbances 
of the stationariness, due to the triggering of aftershock 
sequences, which destroy the universal scaling behavior. 
In order to treat all the distributions in the same way, we 
calculate the rate R only for events in the scaling region, 
i.e., short recurrence times are not considered in the rate. 

The knowledge of the probability distribution of the 
recurrence times allows one to answer two important 
questions about the temporal occurrence of earthquakes. 
First, for a certain region and for M > M c , we can study 
the probability per unit time of an immediate earthquake 
given that there has been a period r without activity, us- 
ing the hazard rate J12[, 
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where the symbol "|" accounts for conditional probabil- 
ity. From the previous formulas we get that A(r) scales 
as A(t) = Rh(Rr), with 
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This function turns out to be monotonically decreasing, 
tending as a power law to the value 1/a as 8 — > oo. So, 
contrary to common belief, the hazard does not increase 
with the elapsed time since the last earthquake, but just 
the opposite; this is precisely the more direct characteri- 
zation of long-term clustering [|l3| . 

Also, one can wonder about the expected time till the 
next earthquake, given that a period To without earth- 
quakes (in the spatial area and range of magnitudes con- 
sidered) has elapsed, 

1 f°° 

e(r = E[t - t | r > t = —— / (r - t )D{t)<1t. 

This function can be referred to as the expected residual 
recurrence time |l2]| and again we find a scaling form for 
it, which is c(tq) = e(Rro) / R, with the scaling function 
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This is an increasing function of 8, which reaches an 
asymptotic value e(9) — > a. Therefore, the residual time 
until the next earthquake should grow with the elapsed 



time since the last one. Notice the counterintuitive be- 
havior that this represents: if we decompose the recur- 
rence time r as t = tq + Tf, with t/ the residual time 
to the next earthquake, the increase of To implies the in- 
crease of the mean value of Tf , but the mean value of r 
is kept fixed. In fact, this is just a more dramatic version 
of the classical waiting-time paradox J14 15|. 



For the particular case of earthquakes this is even more 
paradoxical, since one would say that the longer the time 
one has been waiting for an earthquake, the closer it will 
be, due to the fact that as time passes stress increases on 
the corresponding fault and the next earthquake becomes 
more likely. (Nevertheless, one needs to distinguish be- 
tween earthquakes on a given fault and earthquakes over 
a certain area.) The question was originally put forward 
by Davis et al. fl6|| , who pointed out that if a lognor- 
mal distribution is a priori assumed for the recurrence 
times, the expected residual time increases with the wait- 
ing time. (However, the increase here was associated to 
the update of the distribution parameters as the time 
since the last earthquake, which was taken into account 
in the estimation, increased; we deal with a different con- 
cept of increasing residual time.) Sornette and Knopoff 
0j showed that the increase (or decrease) depends com- 
pletely on the election of the distribution, and studied the 
properties of a number of them. We are going to see that 
the observational data provide direct evidence against 
the simple picture of the next earthquake approaching in 
time. 

Indeed, our mathematical predictions can be con- 
trasted with the catalogs; both the hazard rate and the 
expected residual recurrence time can be directly mea- 
sured with no assumption about their functional form. 
Following their definitions, these quantities are estimated 
as 
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where ?i(ti, T2) denotes the number of quakes with recur- 
rence time in the range (ti,T2) and the sum in c(to) is 
computed only for earthquakes i such that Tj > To (and 
of course M > M c ). From the results displayed in Figs. 1 
(b) and 1 (c) it is apparent that the hazard rate decreases 
with time whereas the expected residual recurrence time 
increases, as we predicted. Moreover, both quantities 
are well approximated by the proposed universal scaling 
functions. 

On the other hand, the part of the recurrence-time dis- 
tribution that accounts for short times displays a typical 
behavior f(8) = Ki/S 1- ", and the corresponding func- 
tions for the survivor function, hazard rate, and expected 
residual return time turn out to be: 
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where the constants Ki and K3 depend on the rest of the 
distribution. An example for these functions with a __ 
0.2 is also represented in Fig. 1, showing the appropriate 
decreasing or increasing tendency in each case. 

For the sake of concreteness, let us consider worldwide 
earthquakes with M > 7.5, which occur at a rate R = 6 
per year, roughly. In the days immediately after one 
event of this type, the expected time to the next one 
(anywhere in the world) is about 2 months (for r = 6 
days, we have Rr = 0.1, and e(0.1) __ 1, see Fig. 1 (c)). 
If after 2 months the quake has not come, the expected 
residual time not only does not decrease but increases to 
2.2 months (e(l) __ 1.1, this would lead to 4.2 months 
between both events), and if the elapsed time reached 
1 year (which is unlikely but not impossible), the ex- 
pect ed waiting time would further increase to 2 , 4 months 
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(e(6 j) ~ 1,2) , — In the same way, the hazard rate would 
drop from 0.7 to 0.5 and to 0.4 month" 1 (h(9) __ 1.4, 1, 
and 0.85), respectively. The same process is reproduced 
at all magnitude and spatial scales in a self-similar man- 
ner. An intuitive explanation of this phenomenon is that 
when the elapsed time since the last earthquake is large, 
the system enters into a long "drought period" in which 
the recurrence time is likely to be very large. Notice 
however that there is no fundamental difference between 
these drought periods and the rest of recurrence times, 
since all of them are governed by the same smooth dis- 
tribution. 

The universality of this behavior demands further ex- 
planation; nevertheless, it suggests the existence of a sim- 
ple mechanism in which, as time passes, the variable that 
triggers rupture runs away from the rupture threshold 
(on average). The "excursions" of this variable would 
keep the memory of the last event stored in the system 
up to very long times to generate the negative aging ob- 
served. 

The considerations reported here should be at the core 
of any research regarding earthquake-occurrence model- 
ing and predictability Jl|,|^,|l|,|| . Finally, in 
order to account for the self-similarity of these processes, 
the concept of self-organized criticality provides the most 
appealing framework up to now |23|j . 
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FIG. 1. Scaling plots of the probability distributions, haz- 
ard rate functions, and expected residual recurrence times for 
all the catalogs analyzed. The values of the occurrence rates 
R are broadly distributed, ranging roughly from 6 year -1 to 1 
hour , so one unit in the horizontal axis can represent from 
1 hour to 2 months. The universal scaling functions fitting 
the data are the ones proposed in the text, with the param- 
eters obtained from the maximum-likelihood estimation; an 
example of fit outside the scaling region is also shown, (a) 
Rescaled survivor functions. The distributions are normal- 
ized for Rt > 0.05, therefore, the left part of the distribu- 
tions does not represent a probability; nevertheless, we have 
considered interesting to show it to illustrate the nonuniversal 
behavior. Times shorter than two minutes are not shown, (b) 
Rescaled hazard rates. The errors at long times are large, due 
to poor statistics, (c) Rescaled expected residual recurrence 
times. Only mean values calculated with 3 or more data are 
displayed. At long times the errors increase even further in 
this e(r) is the difference of two large numbers; nev- 

ertheless, the gaps between the points and the function are 
compatible with the error bars (not shown). 
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